% ECG_FILTER    Generalize function for ECG filters
%
%   [sgofilt,notchfilt,chebyfilt] = ecg_filter(X,Samp_freq)
%   [sgofilt,notchfilt,chebyfilt] = ecg_filter(X)
%
% Composing
% 2. Bandpass filter
% 1. Simple Notch filter
%
% Configuration
% Lowpass:     Chebyshev type II design, adapative orders
% Notch:        Filter out line noise and DC Component
% Smoothing:    Sgolay FIR filter
%
% - Apply Zero-phase-shift filtering
%
% Debug mode:
% When no output is specified, three diagrams will be plotted:
%   1) Frequency Resonse of cascaded IIR filters
%   2) Pole-zero diagram of cascaded IIR filters
%   3) Input and output waveform with spectrograms
%
% See also:     highpass_design,lowpass_design, notch_design, ecg_smooth
% Modified at   27th Oct 2010
% Copyleft by Kelvin,Antony, Karon

function [sgofilt,notchfilt,chebyfilt] = ecg_filter(X,Samp_freq)

if(nargin<2)
     Samp_freq = 553.625;
 end

%% Implement bandpass filter design
%highpass filter
Width=0.01;
[b1 a1] = highpass_design(Width, Samp_freq);

%lowpass filter
lowpass_freq = 80; %corner frequency
Rp = 0.5;
Rs = 40; %attenuation
[b2 a2] = lowpass_design(lowpass_freq, lowpass_freq*1.2,Rp,Rs,Samp_freq);

%%Cascade two filters to make bandpass filter
B_bandpass = conv(b1,b2);
A_bandpass = conv(a1,a2);


%% Implement notch filter design
Width = 0.01;
[b3 a3] = notch_design(50,Width,Samp_freq);

%% Apply filters to the signal
    X_bandpass = filtfilt(B_bandpass,A_bandpass,X);
    chebyfilt = X_bandpass;
    
    X_notch = filtfilt(b3,a3,X_bandpass);
    notchfilt = X_notch;
    
    X_sgo = ecg_smooth(X_bandpass,8,61,Samp_freq);
    sgofilt = X_sgo;

%% Debug mode:
if(~nargout)

 % Cascade Highpass,Lowpass and notch filters
B = conv(B_bandpass,b3);
A = conv(A_bandpass,a3);

figure(1)
%frequency response
freqz(B,A,[],Samp_freq);
title('Frequency Response')

figure(2)
%Pole-zero plot
circle = exp( 1i * linspace(0,2*pi,500) );
plot(roots(B),'bo','markersize',12);
hold on;
plot(roots(A),'rx','markersize',12);
plot(circle,'g');
hold off;
title('Pole-zero Diagram'),xlabel('Real'),ylabel('Imaginary');
grid on;axis equal;

figure(3);
t = (0:length(X)-1)/Samp_freq;

subplot(2,2,1), plot(t,X),xlim([20 23]);
title('Input waveform'),ylabel('Voltage / V');

subplot(2,2,2), specgram(X,[],Samp_freq);
title('Spectrogram'),colorbar;

subplot(2,2,3), plot(t,X_sgo),xlim([20 23]);
title('Output Waveform'),xlabel('time / s'),ylabel('Voltage / V');

subplot(2,2,4), specgram(X_sgo,[],Samp_freq);
title('Spectrogram'),colorbar;

end

